#include "skycalc_real.c"

/* -------------------------------------------------------- */
/* ADDED BY ADRIAN PRICE-WHELAN 2011-06-29 */
/* -------------------------------------------------------- */

double jd2decHour(double jdin) {
	/* given a julian day, prints time only */
	struct date_time date;
	double decHour;
	short dow;
	
	caldat(jdin,&date,&dow);
	decHour = date.h + date.mn/60. + date.s/3600.;
    return decHour;
}

struct almanac {
    double moonIllumination;
    double moonrise;
    double moonset;
    double sunrise;
    double sunset;
    double eveningTwilight;
    double morningTwilight;
    double siderealTimeAtMidnight;
};

struct almanac get_almanac(struct date_time date, 
             double lat, 
             double longit, 
             double elevsea, 
             double elev, 
             double horiz, 
             double stdz, 
             short use_dst) {   
    double jd, jdmid, jdcent, ramoon, decmoon, distmoon;
    double jdb = 0.0;
    double jde = 0.0;
    double sidTimeAtMidnight, sunset, sunrise, moonset, moonrise, eveningTwilight, morningTwilight;
    double moonIllum;
    double geora, geodec, geodist;  /* geocent for moon, not used here.*/
    double locjdb, locjde;
    double rasun, decsun, min_alt, max_alt;
    double hasunset, jdsunset, jdsunrise, sid;
    double hamoonrise, hamoonset, tmoonrise, tmoonset, jdmoonrise, jdmoonset;
    double hatwilight, jdtwilight, jdetw, jdmtw;
    short dow, utOffset; /* day of week */
    float set_to_rise, twi_to_twi, hrs, moonless_hrs, moony_hrs, moon_print;
    struct almanac alm;
    
    // PRINT INPUT PARAMETERS
    //printf("%f\n%f\n%f\n%f\n%f\n%f\n", lat, longit, elevsea, elev, horiz, stdz);
    //exit(0);
    
    find_dst_bounds(date.y, stdz, use_dst, &jdb, &jde);
    locjdb = jdb-stdz/24.;
    locjde = jde-(stdz-1)/24.;
    date.h = 18;  /* local afternoon */
    date.mn = 0;
    date.s = 0;  /* afternoon */

    jd = date_to_jd(date); /* not really jd; local equivalent */
    
    jd = jd + .5;
    jd = jd - 0.25;  /* local midnight */
    jdmid = jd + zone(use_dst,stdz,jd,jdb,jde) / 24.;  /* corresponding ut */
    utOffset = (int) jd2decHour(jdmid);
    
    sidTimeAtMidnight = lst(jdmid,longit);
    //printf("%f", sidTimeAtMidnight);
    
    accumoon(jdmid,lat,sidTimeAtMidnight,elevsea,
             &geora,&geodec,&geodist,&ramoon,&decmoon,&distmoon);    
    lpsun(jdmid,&rasun,&decsun);
    hasunset = ha_alt(decsun,lat,-(0.83+horiz));
    if(hasunset > 900.) {  /* flag for never sets */
        set_to_rise = 0.;
        twi_to_twi = 0.;
        jdcent = -1.;
        goto DO_MOON_ETC;  /* aargh! Ugly flow of control!! But
                   twilight certainly irrelevant if sun up all night, so
           skip it. */
    }
    if(hasunset < -900.) {
        set_to_rise = 24.;
        jdcent = -1.;
        goto CHECK_TWI15;  /* More ugly flow  - 
            checks for twilight even if sun down. */
    }        

    jdsunset = jdmid + adj_time(rasun+hasunset-sidTimeAtMidnight)/24.;  /* initial guess */
    jdsunset = jd_sun_alt(-(0.83+horiz),jdsunset,lat,longit);
    
    //print_time'd this: jdsunset-zone(use_dst,stdz,jdsunset,*jdb,*jde)/24.
    sunset = jd2decHour(jdsunset-zone(use_dst,stdz,jdsunset,jdb,jde)/24.);
    
    jdsunrise = jdmid + adj_time(rasun-hasunset-sidTimeAtMidnight)/24.;
    jdsunrise = jd_sun_alt(-(0.83+horiz),jdsunrise,lat,longit);
    
    //print_time'd this: ((jdsunrise-zone(use_dst,stdz,jdsunrise,*jdb,*jde)/24.),0);
    sunrise = jd2decHour(jdsunrise-zone(use_dst,stdz,jdsunrise,jdb,jde)/24.);
    
    if((jdsunrise > 0.) && (jdsunset > 0.)) {
        set_to_rise = (jdsunrise - jdsunset) * 24.;
        jdcent = (jdsunrise + jdsunset) / 2.;
    }
    else {
        oprntf("Sunrise not correctly computed.");
        jdcent = -1.;
    }

    /* This block checks for and prints out 15-degree twilight
       as appropriate. */

    CHECK_TWI15: hatwilight = ha_alt(decsun,lat,-15.);
    if(hatwilight < -900.) {
        twi_to_twi = 24.;
        goto DO_MOON_ETC;  /* certainly no 12-degree twilight */
    }
    if(hatwilight > 900.) {
        twi_to_twi = 0.;
    }               

    /* compute & print evening twilight */
    jdtwilight = jdmid + adj_time(rasun+hatwilight-sidTimeAtMidnight)/24.;  /* rough */
    jdtwilight = jd_sun_alt(-15.,jdtwilight,lat,longit);  /* accurate */
    jdetw = jdtwilight;
    if(jdtwilight > 0.) {
        // Evening twilight: print_time((jdtwilight-zone(use_dst,stdz,jdtwilight,*jdb,*jde)/24.),0);
        eveningTwilight = jd2decHour(jdtwilight-zone(use_dst,stdz,jdtwilight,jdb,jde)/24.);
        
        /*sid = lst(jdtwilight,longit);
        print_tz(jdtwilight,use_dst,*jdb,*jde,zabr);
        oprntf(";  LMST at evening twilight: ");
        put_coords(sid,0,0);
        oprntf("\n");*/
    }

    /* Now do morning twilight */
    jdtwilight = jdmid + adj_time(rasun-hatwilight-sidTimeAtMidnight)/24.;
    jdtwilight = jd_sun_alt(-15.,jdtwilight,lat,longit);
    jdmtw = jdtwilight;
    if(jdtwilight > 0.) {
        // Morning twilight: print_time((jdtwilight-zone(use_dst,stdz,jdtwilight,*jdb,*jde)/24.),0);
        morningTwilight = jd2decHour(jdtwilight-zone(use_dst,stdz,jdtwilight,jdb,jde)/24.);
        
        /*sid = lst(jdtwilight,longit);
        print_tz(jdtwilight,use_dst,*jdb,*jde,zabr);
        oprntf(";  LMST at morning twilight: ");
        put_coords(sid,0,0);*/
    }

    if((jdetw > 0.) && (jdmtw > 0.)) twi_to_twi = 24. * (jdmtw - jdetw); 


    DO_MOON_ETC: 
    /* NIGHT CENTER
    if(jdcent > 0.) {
        oprntf("night center: ");
        print_time((jdcent-zone(use_dst,stdz,jdcent,*jdb,*jde)/24.),0);         
        print_tz(jdcent,use_dst,*jdb,*jde,zabr);
    }
    oprntf("\n\n");*/
    
    min_max_alt(lat,decmoon,&min_alt,&max_alt);  /* rough check -- occurs? */       
    if(max_alt < -(0.83+horiz)) {
        //oprntf("Moon's midnight position does not rise.\n");
        jdmoonrise = -1.;
        goto MORE_MOON;
    }
    if(min_alt > -(0.83+horiz)) {
        //oprntf("Moon's midnight position does not set.\n");
        jdmoonrise = 1.;
        goto MORE_MOON;
    }                       

    /* compute moonrise and set if they're likely to occur */

    hamoonset = ha_alt(decmoon,lat,-(0.83+horiz)); /* rough approx. */
    tmoonrise = adj_time(ramoon-hamoonset-sidTimeAtMidnight);
    tmoonset = adj_time(ramoon+hamoonset-sidTimeAtMidnight);
    jdmoonrise = jdmid + tmoonrise / 24.;
    jdmoonrise = jd_moon_alt(-(0.83+horiz),jdmoonrise,lat,longit,elevsea);
    jdmoonset = jdmid + tmoonset / 24.;
    jdmoonset = jd_moon_alt(-(0.83+horiz),jdmoonset,lat,longit,elevsea);
    
    if(fabs(set_to_rise) > 10.) moon_print = 0.65*set_to_rise;
    else moon_print = 6.5;  

    /* more lousy flow of control, but it's nice to see the event which happens
       first printed first .... */

    if((jdmoonrise > 0.) && (fabs(tmoonrise) < moon_print)) { 
        // MOON RISE print_time((jdmoonrise-zone(use_dst,stdz,jdmoonrise,*jdb,*jde)/24.),0);
        moonrise = jd2decHour(jdmoonrise-zone(use_dst,stdz,jdmoonrise,jdb,jde)/24.);
    }

    if((jdmoonset > 0.) && (fabs(tmoonset) < moon_print)) {
        // MOON SET print_time((jdmoonset-zone(use_dst,stdz,jdmoonset,*jdb,*jde)/24.),0);
        // print_tz(jdmoonset,use_dst,*jdb,*jde,zabr);
        moonset = jd2decHour(jdmoonset-zone(use_dst,stdz,jdmoonset,jdb,jde)/24.);
    }
    goto MORE_MOON;

    MORE_MOON: moonIllum=0.5*(1.-cos(subtend(ramoon,decmoon,rasun,decsun)));
    
    alm.moonIllumination = moonIllum;
    alm.moonrise = fmod(moonrise + 7.0, 24.0);
    alm.moonset = fmod(moonset + 7.0, 24.0);
    alm.sunrise = fmod(sunrise + 7.0, 24.0);
    alm.sunset = fmod(sunset + 7.0, 24.0);
    alm.eveningTwilight = fmod(eveningTwilight + 7.0, 24.0);
    alm.morningTwilight = fmod(morningTwilight + 7.0, 24.0);
    alm.siderealTimeAtMidnight = sidTimeAtMidnight;
    
    return alm;
}


